function [dfile, data, err] = read_acq_mat(filename, RL)
%
% read matlab acquisition file structures. Returns data for
% records in RL list from the file.
% and fills the dfile structure.
% Operation is exact, but not necessarily inspired.
%
% 10/15/2000
% Paul B. Manis, Ph.D.
% pmanis@med.unc.edu
%
% some speed enhancemenets 2/16/05 P. Manis.
% using newer matlab constructs, as in dynamic field names.
%
% Updates to fix multichannel reading and clean up code 1/13/2009
% Paul B. Manis

err = 1;
data = [];
if(nargin  < 2 || nargin > 2)
    fprintf(1, 'read_acq_mat: requires exactly 2 arguments\n');
    return;
end;


[path, fname, ext] = fileparts(filename);
dfile = init_dfile;
dfile.fullfile = filename;
dfile.filename = fname;
dfile.path = path;
dfile.ext = ext;

hdr = load(filename, 'HFILE'); % read the header to get top information
if(isempty(hdr) || isempty(fieldnames(hdr)))
    fprintf(1, 'read_acq_mat: File %s lacks Header information\n', filename);
    return;
end;
if(~strmatch(fieldnames(hdr), 'HFILE'))
    fsprintf(1, 'read_acq_mat: HFILE in hdr not found in %s\n', filename);
    return;
end;
hdr = hdr.HFILE;
comment = char(hdr.Experiment.v);
dfile.comment = [];
for j = 1:size(comment,1)
    dfile.comment=strcat(dfile.comment, sprintf(' %s', comment(j,:)));
end;
dfile.Signature = hdr.Signature.v;
dfile.mode = 5; % mat file mode - not used

warning off
try
    index = load(filename, 'Index');
catch
    if(isempty(RL))
        fprintf(1, 'read_acq_mat: ***%s has empty Index\n', filename);
    else
        fprintf(1, 'read_acq_mat: Index is empty in %s', filename);
    end;
    return;
end;
warning on;

if(isempty(index) || ~isstruct(index) || isempty(fieldnames(index)))
    if(isempty(RL))
        fprintf(1, 'read_acq_mat: %s has empty Index\n', filename);
    else
        QueMessage(sprintf('read_acq_mat: Index is empty in %s', filename), 1);
    end;
    return;
end;
while(isempty(strmatch('date', fieldnames(index), 'exact')))
    index = index.(char(fieldnames(index)));
end;

if(~isempty(RL))
    dfile.record = min([index.record]); %should be 1 for first record
    dfile.nr_channel = 2; % use the default; it can change anytime
end
dfile.records_in_file = max([index.record]);

if(isempty(RL)) % catch when we just read the header for information
    err = 0;
    return;
end
if (max(RL) > dfile.records_in_file)
    QueMessage(sprintf('read_acq_mat: Last record (%i) greater than length of file (%i recs in file)',max(RL),dfile.records_in_file), 1);
    return
end

data_blks = strmatch('DATA', char(index.type), 'exact');
if(isempty(data_blks))
    QueMessage('read_acq_mat: No data blocks in index? ', 1);
    return;
end;

% make a fake end block....
maxind = length(index)+1;
index(maxind).record=10000000;
data_blks = [data_blks; maxind];

dfile.frec=RL(1);
dfile.lrec=RL(end);
dfile.record = [];
dfile.rate = [];
dfile.gain = [];
dfile.rate = [];
dfile.nr_points = [];
dfile.block = data_blks; % save the data blocks too (only in this version)
dfile.dblock = [];
dfile.ztime = [];
rec_ptr = 1;
S = load(filename, '-regexp', '[sd]f[0-9]');
for indxi = 1:length(data_blks) % look for the data in all the blocks (it doesn't hurt to search)
    ind = data_blks(indxi);
    if(any(RL >= index(ind-1).record) && any(RL <= index(ind).record - 1)) % any records in this block?
        this_block = index(ind).block;
        if(~isempty(this_block))
            sfc = sprintf('sf%d', this_block);
            dfc = sprintf('df%d', this_block);
            sfile =S.(sfc);
            df=S.(dfc); % using new dynamic field name method.
        else
        end
        dfile.dblock = [dfile.dblock this_block]; % save the blocks read
%
% build channel list order from df.ChannelList.v (new format)
% 
        ch_order = [1:length(df.Channels.v)]; % default order is 1 to n
        if isfield(df, 'ChannelList') 
            list = df.ChannelList.v; % get the list
            switch( upper(char(df.Data_Mode.v)) )
                case 'VC'
                case 'CC'
                    if list(1) == 'i'
                        o1 = ch_order(1);
                        ch_order(1) = ch_order(2);
                        ch_order(2) = o1;
                    end;
                otherwise
            end;
            
        end   
        QueMessage(sprintf('Reading Matfile data, block = %d', this_block));
        nm = strmatch(sprintf('db_%d', this_block), char(index.MatName), 'exact'); % only one can match.
        v=load(filename, index(nm).MatName); % load all the relevant data
        if(isempty(v))
            QueMessage(sprintf('read_acq_mat.m: There appears to be no data in block db_%d', this_block));
            return;
        end;

        fn=fieldnames(v); % find out the names in the list
        d = strmatch('db_', fn); % search for data blocks
        if(isempty(d))
            QueMessage(sprintf('read_acq_mat.m: There appears to be no data in block db_%d', this_block));
            return;
        end;
        d2 = v.(fn{d});
        if(isempty(d2))
            QueMessage(sprintf('read_acq_mat.m: Data in block %d did not evaluate', this_block));
            return;
        end;
        drl = get_record_list(d2); % find the record numbers present in the block
        BRL = intersect(RL, drl); % access the records matching those we request
        if(~isempty(BRL))
            dfile.rate = [dfile.rate df.Sample_Rate.v*ones(1, length(BRL))];
            dfile.nr_points = df.Points.v;
            if(indxi == 1) % preallocate the data area itself!
                data = zeros(length(RL), dfile.nr_points);
            end;
            dfile.ftime = 0;
            switch(upper(char(df.Data_Mode.v)))
                case 'VC' % put the for loops inside the switch for speed.
                    dfile.dmode = 'VC';
                    for i = 1:length(BRL)
                        indx = BRL(i)-drl(1)+1;
                        dn = fieldnames(d2{indx});
                        dlen = dfile.nr_points;
                        adat=d2{indx};
                        pd = (1:dlen);
                        for nc = 1:length(df.Channels.v)
                            data(rec_ptr,pd+(nc-1)*dlen) = adat.(dn{1}).data(:,nc);
                        end;    
                        ev = d2{indx}.(dn{1}).events; 
                        k = strmatch('Trigger', {ev.Type});
                        if(~isempty(k) && k > 0)
                            ti = ev(k).Data.AbsTime;
                            dfile.ztime(rec_ptr) = ti(4)*3600+ti(5)*60+ti(6); % express time in seconds.
                        else
                            dfile.ztime(rec_ptr)=i;
                        end;
                        dfile.ftime(rec_ptr) = 0; % not used in new structure
                        dfile.record(rec_ptr) = BRL(i);
                        dfile.channels(rec_ptr) = length(df.Channels.v);
                        dfile.nr_channel(rec_ptr) = dfile.channels(rec_ptr);
                        dfile.rate(rec_ptr)=df.Sample_Rate.v;
                        dfile.gain(rec_ptr,:)=ones(1,8);
                        dfile.low_pass(rec_ptr,1:length(df.Low_Pass.v))=df.Low_Pass.v;
                        rec_ptr = rec_ptr + 1;
                    end;
                case 'CC' 
                    dfile.dmode = 'CC';
                    getfromsfile = 0;
                    if(getfromsfile)
                        [pp] = fileparts(which('datac'));
                        h = fopen(slash4OS([append_backslash(pp) 'temp.m']), 'w'); % save the method so we can calculate on the fly
                        mc = sfile.method_code;
                        fprintf(h, '%s', mc');
                        fclose(h);
                    end;

                    for i = 1:length(BRL)
                        indx = BRL(i)-drl(1)+1;
                        adat=d2{indx};
                        % convert current time to minutes with decimal seconds
                        % and store in ftime.
                        dfile.record(rec_ptr) = BRL(i);
                        dfile.channels(rec_ptr) = length(df.Channels.v);
                        dfile.nr_channel(rec_ptr) = dfile.channels(rec_ptr);
                        dfile.rate(rec_ptr)=df.Sample_Rate.v;
                        dfile.gain(rec_ptr,:)=ones(1,8);
                        dfile.low_pass(rec_ptr,1:length(df.Low_Pass.v))=df.Low_Pass.v;

                        dn = fieldnames(d2{indx});
                        dlen = dfile.nr_points;
                        pd = (1:dlen);
                        data(rec_ptr,pd) = adat.(dn{1}).data(:,ch_order(1));
                        switch(getfromsfile)
                            case 0
                                data(rec_ptr,pd+dlen) = adat.(dn{1}).data(:,ch_order(2));
                            case 1 % for the axopatch amplifiers, the intended command is a better choice.
                                p=df.Sample_Rate.v*dfile.nr_channel(i)/1000;
                                tbx=0:p:(dfile.nr_points-1)*p;
                                QueMessage( 'Computing Stimulus waveform...', 1);
                                eval('[istim stb out_rate err   ] = temp(sfile);');
                                data(rec_ptr,dlen+1:2*dlen) = interp1(double(stb{i}.v), double(istim{i}.v), tbx); %#ok<USENS>
                                dfile.stim{i} = data(rec_ptr, dlen+1:2*dlen);
                            case 2
                                p=df.Sample_Rate.v*dfile.nr_channel(i)/1000;
                                tbx=0:p:(dfile.nr_points-1)*p;
                                QueMessage('Getting stimulus waveform from data', 1);
                                data(rec_ptr,dlen+1:2*dlen) = interp1(double(stb{1}.v), double(istim{i}.v), tbx);
                                dfile.stim{i} = data(rec_ptr, dlen+1:2*dlen);
                        end;
                        if(dfile.channels(rec_ptr) > 2)
                            for ich = 3:dfile.channels(rec_ptr)
                                data(rec_ptr,pd+(ich-1)*dlen) = adat.(dn{1}).data(:,ch_order(ich));
                            end;
                        end;

                        ev = d2{indx}.(dn{1}).events; 
                        k = strmatch('Trigger', {ev.Type});
                        if(k > 0 && ~isempty(k))
                            ti = ev(k).Data.AbsTime;
                            dfile.ztime(rec_ptr) = ti(4)*3600+ti(5)*60+ti(6); % express time in seconds.
                        else
                            dfile.ztime(rec_ptr)=i;
                        end;
                        rec_ptr = rec_ptr + 1;
                    end;
                otherwise
            end;
        else
            QueMessage(sprintf('read_acq_mat: Record list for block %d empty. Possible File Corruption!', this_block), 1);
        end;
    end;
end;

if(isempty(dfile.rate))
    QueMessage(sprintf('read_acq_mat: No data loaded? Probable File Corruption'));
    return;
end;

% now fill the header information
dfile.slow=0; % not used in current version
dfile.durs = [5 100 0 0];
if(exist('sfile', 'var') && strcmp(sfile.Method.v, 'steps'))
    dfile.durs = sfile.Duration.v; % also store the durations from the last stim...
    dfile.seq = sfile.Sequence.v;
    dfile.seqpar = sfile.SeqParList.v;
    dfile.seqstep = sfile.SeqStepList.v;
    dfile.levs = sfile.Level.v;
    w = sfile.waveform;
    fnw = fieldnames(w{1});
    if(strmatch(fnw, 'v2sco'))
        shkl=find(w{1}.v2sco ~= 0);
    else
        shkl = find(w{1}.v2 ~= 0);
    end;
    if(~isempty(shkl))
        dfile.shkl = (shkl(1)+1)*1000/sfile.outrate;
    else
        dfile.shkl = []; % no detected shock latency.
    end;

end;

if(~isempty(RL))  % we sometimes override these with the ctl structure...
    dfile.refvgain=1;
    dfile.refigain=1;
    dfile.refwgain=1;
end;
%profile viewer
err = 0;
return


function [rl] = get_record_list(d)
rl = zeros(length(d), 1);
for i = 1:length(d)
    dn = char(fieldnames(d{i})); % get the fieldname of the entry
    j = find(dn == '_'); % parse data to find record name
    sn = dn(j(end)+1:end); % get the record
    rl(i) = str2double(sn); % make a number
end;
return;
